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ABSTRACT 

Context. Theory predicts that low-mass protoplanets in a protostellar disc migrate into the central star on a time scale that is short 
compared with the disc lifetime or the giant planet formation time scale. Protoplanet eccentricities of e ^ H/r can slow or reverse 
migration, but previous 2D studies of multiple protoplanets embedded in a protoplanetary disc have shown that gravitational scattering 
cannot maintain significant planet eccentricities against disc-induced damping. The eventual fate of these systems was migration into 
the central star. 

Aims. Here we simulate the evolution of low-mass protoplanetary swarms in three dimensions. The aim is to examine both protoplanet 
survival rates and the dynamical structure of the resulting planetary systems, and to compare them with 2D simulations. 
Methods. We present results from a 3D hydrodynamic simulation of eight protoplanets embedded in a protoplanetary disc. We also 
present a suite of simulations performed using an A'-body code, modified to include prescriptions for planetary migration and for 
eccentricity and inclination damping. These prescriptions were obtained by fitting analytic formulae to hydrodynamic simulations of 
planets embedded in discs with initially eccentric and/or inclined orbits. 

Results. As was found in two dimensions, differential migration produces groups of protoplanets in stable, multiple mean-motion 
resonances that migrate in lockstep, preventing prolonged periods of gravitational scattering. In almost all simulations, this leads to 
large-scale migration of the protoplanet swarm into the central star in the absence of a viable stopping mechanism. The evolution 
involves mutual collisions, occasional instances of large-scale scattering, and the frequent formation of the long-lived, co-orbital 
planet systems that arise in > 30% of all runs. 

Conclusions. Disc-induced damping overwhelms eccentricity and inclination growth due to planet-planet interactions, leading to 
large-scale migration of protoplanet swarms. Co-orbital planets are a natural outcome of dynamical relaxation in a strongly dissipative 
environment, and if observed in nature would imply that such a period of evolution commonly arises during planetary formation. 

Key words, planet formation- extrasolar planets- -orbital migration-protoplanetary disks 

1. Introduction One of several problems associated with this picture is the 

rapid inward migration experienced by a protoplanet due to the 

O The observed lower limits on extrasolar planet mass are contin- gravi t ational interac tion between it and the gaseous disc jWar^ 

uing to decrease, while at the same time the number of known I1997I: iTanakaet al.ll2002h . In particular, understanding the for- 

--j multi ple planet s ystgins is continuing to grow (Rivera et al] nation of giant planets, which must spend more that lO*" years 

X; ; |2005t iLovis et al.||2006|; |Udry et al.||2007D. Multiple planet sys- in the disc in order to accrete a gas envelope, remains an un- 

H ; tems containing sub-Neptune mass planets are now being dis- solved problem. Their sohd cores have migration times shorter 

- - . covered. Missions such as CoRoT and Kepler are expected to than both the gas accretion time scale and the disc life time 

discover further sub-Neptune mass planets beyond the reach of (Pa paloizou & NelsonlllOosh . Referred to as type I migration, 

current observations, leading the way toward finding planetary this drift makes it hard to understand how gas giants can form 

systems more hke our own. at all without being accreted by the central star. Solving the type 

Formation of a planetary system is beheved to involve ac- I problem remains an active area of research, and recently sug- 

cretion within a protoplanetary disc, involving essentially three gested remedies incliide: stoc hastic m igration in a turbulent disc 

steps: coagulation of dust grains into small (~ 1 km) planetes- dNelson & Papaloizoul l2004t iNelsonI 12005,) ; corotation torques 

imals; runaway growth of planetesimals into larger (~ 100- o perating in a region of positive gradient in disc surface density 

1000 km) protoplanets; oligarchic growth by planetesimal ac- (Masset et al. 2 00^; corotation torque s in radiatively inefficient 

cretion into larger planetary cores. Those cores forming beyond discs (jPaardekooper & Mellemal l2006'). In this paper we further 

the snow line are expected to reach masses of ~ 10 M®, accreting explor e the role of protoplanet scattering tCresswell & NelsonI 

a gaseous envelope to become gas giants if they form before disc ( l2006h — hereafter Paper I), 

dispersal, or ice giants should they form late in the disc's lifetime ,, . 

jBodenheimer & PoUackl [1981 IPoUack et alJ [19961) . Smaller, , ^^^^ ^^^^^^"^ have examined the interactions of multiple 

T,/r , f ,. , • »u ■ *u * planetary embryos, or fully formed planets, within protoplane- 

Mars-mass bodies result from oligarchic growth in the terres- ^ ,. xr • ,• .• , . 

.., , • tary discs. Numerical simulations of the oligarchic growth phase 

trial zone, and accrete via giant impacts to form an inner system , ... ■ r , 

f , ^ ^ • , , ^ jj^r — r — p -iir — Tni inno k have been used to examine the interactions of planetary em- 

of rocky, terrestrial-mass planets ^Chambers & Wetherill 1998ft . , jrr . , n ,. ^^^^ ^, , x, . ; . 1 

^ ^ ^ ^' ^ bryos dKpkubo & Ida 2000; Thommes et al. 2003). McNeiretap 

(I2OO5I) found that a proto-terrestrial system may form against 

Send offprint requests to: RCresswell@qmul.ac.uk type I migration by enhancing the disc mass by a factor of 2-4. 
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Resonant captur e between two planets in the 1 -20 range has 
been s tudied by 'Papaloizo u & Szuszkiewicj ( l2005h . iThommesI 
( l2005b suggests that the formation of the first giant planet may 
be a significant step in the formation of a planetary system, by 
capturing smaller cores in resonance and preventing further type 
I migration. Previously, pairs of giant planets in resonance have 
been examined, often with direct application to a specific system 
such as GJ 876 (Snellgrove et al. 2001; Kle y et al. 2005). 

One area that has not been addressed to the same depth is 
the issue of how a swarm of protoplanetary cores, of Earth mass 
and above, will evolve under the influence of a surrounding pro^ 
toplanetary disc. Models of oligarchic growth ( K okubo & Idal 
l2000t iThommes et al.l 12003 ) predict that a number of cores 
should form coevally, separated in radius by ~ 8 mutual Hill 
radii. Diff'erential type I migration may cause these bodies to un- 
dergo close encounters, leading to gravitational scattering and 
the pumping of eccentricitie s. 

iPapaloizou & LarwoodI (l2000h found that type I migration 
can be slowed or reversed when a protoplanet embedded in a 
disc achieves an eccentricity e <: l.lH/r, where H/r is the disc's 
scale height-to-radius ratio, raising the possibility of the mutual 
interactions of a swarm of cores sustaining significant eccentric- 
ities and slowing type I migration for at least some of them. In 
an earlier study, however, using 2D A^-body and hydrodynamic 
models we found that the damping action of the disc is too strong 
to sustain eccentricities in this manner (see Paper I). Instead the 
cores form resonant groups and migrate inwards in lockstep. Due 
to the 2D modelling it is possible that the collision probability 
was overestimated, removing a higher number of protoplanets 
than should be properly expected, and the influence of planetary 
inclinations on reducing the migration rate was neglected. 

To address these issues, we perform 3D numerical simula- 
tions of swarms of planetary cores embedded in a protoplanetary 
disc. In the first instance we present the results of a full, 3D hy- 
drodynamic simulation of a protoplanetary disc containing eight 
protoplanets. This simulation provides a good illustration of the 
early stages of evolution, but cannot be run for long times. To ex- 
amine the long term evolution, we also present 3D simulations 
performed with an A^-body code, modified to include prescrip- 
tions for migration, and eccentricity and inclination damping. 
These prescriptions were obtained by fitting analytic formulae 
to numerous 3D hydrodynamic simulations of protoplanets on 
eccentric and/or inclined orbits embedded in a protoplanetary 
disc. The results of these 3D multiplanet simulations suggest that 
gravitational interactions among a swarm remains ineff'ective at 
maintaining the requisite eccentricities to slow/stop migration, 
due to the strong damping from the disc. We also find that many 
co-orbital planets form naturally from such a migrating popula- 
tion, raising the possibility of their detection among the observed 
extrasolar hot Neptunes and super-earths. 

The plan of this paper is as follows. In Sect. 2 we describe 
the basic equations of motion. In Sect. 3 we describe the two 
numerical schemes and provide equations describing the disc's 
action on the protoplanets. In Sect. 4 we describe the initial con- 
ditions. In Sect. 5 we present a hydrodynamic multiple-planet 
model. In Sect. 6 we present the results obtained from the mod- 
ified A^-body scheme, and discuss the trends and implications of 
these results. We give our conclusions in Sect. 7. 



2. Equations of Motion 

An unperturbed protoplanetary disc with constant aspect ratio 
H/r can be conveniently described using spherical polar coor- 



dinates (r,9,^) with the origin located at the central star The 
continuity equation is expressed as 



^ + V ■ (pv) = 0, 
ot 



(1) 



and the three components of the momentum equation are written 



dipVr) , ^ , , P(^9 + dp 

+ V ■ (PV,V) = — - P— + fy (2) 

ot r or Or 



dipvg) pvyve Pi^VCOtfl 

+ V ■ (pVgV) = -I- 

at r r 



I dp pd<b 
r dr r do 



(3) 



PVVV0 pVgV^ cote 



dp 



p 



50 



r sin 6 d(f) r sin i 



(4) 



Here p denotes the density, p is the gas pressure and f^, fg, 
are the viscous forces per unit volume in the radial, meridional 
and azimuthal directions, respectively. Vr, vg and denote the 
corresponding velocities. The gravitational potential O is given 
by 




(5) 



where M, is the stellar mass, e is a softening parameter, and the 
summations are over all protoplanets p with masses ntp. The sub- 
script '/?' denotes evaluation at the location of the protoplanet. 
The latter two terms in Eq. |5]result from acceleration of the co- 
ordinate system due to the gravity of the protoplanets and the 
protostellar disc. The integral is performed over the volume of 
the disc. 

Each protoplanet experiences the gravitational acceleration 
from the central star, the other protoplanets and the protostellar 
disc. The equation of motion for each protoplanet is; 



dVp 
dt2 

where 



-VOd 



-Z 



Gnipi 



z 



Gnin' 



= -G 



p(r)dr 



+ G 



X 



dm(r) 



(6) 



(7) 



is the gravitational potential due to the disc. The integrals are 
again performed over the volume of the disc. The final terms in 
Eqs.|6]and|2]are the indirect terms arising from the acceleration 
of the coordinate system by the protoplanets and disc, respec- 
tively. 

We adopt a locally isothermal equation of state such that the 
pressure and density are related by 



(8) 
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where Cs - (7)^k is the isothermal sound speed, a fixed func- 
tion of distance from the central star vk is the local Keplerian 
velocity. All disc models have aspect ratio H/r - 0.05. We also 
adopt the alpha model for the disc viscosity such that the kine- 
matic viscosity v - acsH, with a = 5 x 10"^. 

3. Numerical methods 

We use two distinct numerical schemes: a 3D hydrodynamic disc 
model together with embedded planets is computed using a hy- 
drodynamics code (NIRVANA); a suite of simulations are com- 
puted using a much faster A^-body code which has been adapted 
to emulate the effects of orbital migration, and eccentricity and 
inclination damping on the protoplanets due to the protoplane- 
tary disc. We describe each in turn in the following sections, and 
demonstrate that the modified A^-body code agrees well with re- 
sults from the hydrodynamic code. 

3.1. Hydrodynamic scheme 

We use a modified version of the grid-based code NIRVANA 
to co nduct 3D hydrodynamical simulations ( Zeigler & Yorke, 
Il997h . The code has previously been applied to a variety of 
disc-planet numerical studies in both two and three dimensions. 
Furthe r details may be found in iNelson et al.l (l2000h . ICresswell 
(120061) . ICresswell et al.l(l2007l) and in Paper I. 

The motion of the protoplanets is int egrated using a 5*-order 
Runge-Kutta scheme dPress et al As in Paper I disc self- 

gravity is neglected. We employ reflecting boundary conditions 
in the meridional direction, and wave-damping conditions at the 
inner and outer radial boundaries. The implementation of these 
is described in Paper I, along with a description of the time step 
control procedure. 

In commo n with many hydrodynamical disc simulations and 
ICresswell et al. (2007), we adopt a scale height Hjr - 0.05 and 
surface density profile S oc r""^. A disc opening half-angle of 
10° then models three and a half scale heights. Our adoption of a 
locally isothermal equation of state means that wave propagation 
is primarily confined to the radial direction, such that the use of 
reflecting boundary conditions at the meridional boundaries does 
not lead to significant wave reflection toward the disc midplane. 
The development of inclined orbits for the protoplanets, how- 
ever, can cause the excitation of bending waves in the disc, which 
could in principle be affected by the vertical boundary condi- 
tions. Test simulations of inclined planets in discs using zero- 
gradient outflow boundary conditions, hwoever, indicate that the 
results presented in this paper are not strongly affected by the 
choice of meridional boundary conditions. The adoption of re- 
flecting conditions prevents the slow loss of mass from the disc 
which accompanies the use of open boundaries. 

3.1.1. Code resolution and modelling 

Due to the large computational expense of 3D disc-planet stud- 
ies, the spatial resolution used in each model must be care- 
fully weighedjigainst the durat;ion and number of simulations 
required. In lCresswell et al.l (l2007h it was found that the coarsest 
resolution tested produced results of a similar character in mi- 
gration rate and eccentricity and inclination damping rates to the 
highest resolutions tested, with some difference in absolute mi- 
gration/damping times; numerical convergence of the code was 
established with these high-resolution models. Since we have 
needed to run many simulations to obtain fitting formulae for the 



modified A^-body code, we have chosen a resolution for those 
simulations that leads to accurate results, but which makes the 
problem tractable: (Nr,Ng,N^) = (128,40, 300). This resolution 
corresponds closely t o the lower-resolution runs presented in 
ICresswell et al.l (|2007|) . The limits of the computational domain 
are defined by flie intervals r = [0.4, 2.5], = [80°, 100°], = 
[0, 27r]. With the set-up described, we follow the evolution of 
a 10 Earth mass (M®) planet on a variety of eccentric and in- 
clined orbits, and fit the migration and eccentricity and inclina- 
tion damping rates using simple nonlinear functions of e and /. 

We also performed one hydrodynamic simulation consist- 
ing of eight protoplanets embedded in a protoplanetary disc. 
The resolution adopted for this simulation was (Nr,Ng,N^) = 
(288, 40, 444), with the limits of the computational domain given 
by r = [0.6, 3.0], 61 = [80°, 100°], ^ = [0, 2n]. 

3.2. N-body scheme 

The second method we use to evolve a disc-planet system is an 
A^-body integrator (to which we apply the name HENC-3D) with 
analytic functions which model the effects of type I migration, 
and eccentricity and inclination damping. The integrator used is 
the same 5*-order Runge-Kutta routine as used in NIRVANA. 
To prevent the time step from becoming too low, protoplanets 
were removed from the simulation if they fell within 10 Solar 
radii of the origin. 

Due to the computational expense of 3D hydrodynamic com- 
putations, this is our primary method of evolving multiple-planet 
systems; even if it were practical to model a disc in NIRVANA 
over the scales required, HENC-3D is ~ 10^ times faster at 
evolving such systems even at the moderate resolution we em- 
ploy. 

3.2.1. Prescriptions for eccentricity and inclination damping 
and migration 

We construct prescriptions for the migration rates and eccentric- 
ity and inclination damping rates which are incorporated into 
the A^-body code. These allow us to follow the evolution of pro- 
toplanet swarms for longer than is possible with the hydrody- 
namic code. Our starting point is the formulae for damping rates 
obtained by Tanaka & Ward (2004) and the migration rates ob- 
tained bv ITanaka et ah (2002), which we modify using multi- 
plicative factors that account for changes in damping and migra- 
tion rates due to planetary eccentricity and/or inclination. These 
multiplicative factors are obtained by fitting formulae to the re- 
sults of numerous hydrodynamical simulations, performed with 
different values of e and /. The results of a subset of these sim- 
ulations are shown in Figs. 1-3. In the case of eccentricity and 
inclination damping, the formulae are found ed on the damping 
time scale derived bv lTanaka & WardI ( l2004i) . and given by their 
Eq. 49: 

fwave = — =^ -U^p' (9) 

where Q. - v^jr is the angular velocity of the unperturbed disc 
and S is the surface density 



•J —cx 



p dz. 



(10) 



The eccentricity damping time that best fits our simulation re- 
sults is given by: 

fwave ^^^^ 



(■wave 



0.780 
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1-0.141 — 

\Hlr 



+ 0.06 



HIr 



+ 0.18 



H/rj\H/r 



where the terms in square brackets provides the modification 
to the damping rate caused by the finite eccentricity and/o r in- 
clination. As was fou nd by iPapaloizou & LarwoodI (l2000h and 
ICresswell et al.l (l2007h . our results show that de/dt cc e'^ at high 
e, with exponential decay for e < Hjr. The inchnation damping 
time, fi, is given by: 



_ 'wave 

^' " 0.544 

1 - 0.30 



(12) 



, +0.24(— ] +0.14(— ] , 
HIrj \H/rj \H/rj \H/r 



where again the terms in square brackets modify the damp- 
ing rate b ecause of finite e ccentricity and/or inclination. As de- 
scribed in ICresswell et al.l ( |2007 ) we find di/dt oc at high /, 
with exponential decay for ; < H/ r. 

The prescript ion for the ni igration time, fm, is based on the 
results of Tanaka et a D (120021) . who considered circular orbits 
only. The results of our simulations are well-fitted by the the 
expression: 



2f„ 



2.7+ l.lyS 



ive / H\ ^ 

\.\B\ r I 



where 



!0.070|— +0.085 — 
\H/rj \H/r 



^ (2.25///r) {iMH/r) 



\P{e)\ 

\4 



(13) 



• 0.080 



H/rl\H/r 



1 



\2.02H/r) 



and jS is given by !E(r) oc r"^. The sign dependency on P(e) 
allows torque reversal at s ufficiently high e. We note that 
IPapaloizou & LarwoodI (1200 0) undertook a study of the disc- 
induced migration and eccentricity damping of low-mass pro- 
toplanets using linear theory, and derived similar expressions to 
those given in Eqs.[TT]and[T3] Our expressions differ from theirs 
only because we found that Eqs. (TTj and [T3] better fit our hy- 
drodynamical simulations, and we consider both eccentric and 
inclined orbits. 

We implement the following expressions in the A^-body code 
as accelerations experienced by the protoplanets due to the disc, 
using values of t^, and fi obtained from Eqs. nTHlJ] 



"m = , 

fm 



(v.r)r 



(14) 

(15) 
(16) 



where k is the unit vector in the z-direction. 



3.2.2. Comparison between A^-body and hydrodynamic code 

We briefly demonstrate the agreement between the modified A^- 
body and hydrodynamic codes. Figure [1] shows the orbital evo- 
lution of a planet on a variety of eccentric orbits, and the corre- 
sponding planetary migration. The decay of eccentricity at both 



NIRVANA data 
HENC-3D model _ 





Fig. 1. {Top) Eccentricity evolution of a 10 Me protoplanet on 
a variety of eccentric orbits. Solid lines are the results from 
NIRVANA, crosses from the A^-body code. {Bottom) The cor- 
responding migration rates. 



high and low eccentricity is modelled accurately, with the largest 
deviation at the transition from quadratic to exponential decay. 
We find that the exchange of angular momentum between disc 
and planet reverses sig n when e ^ 2j//r, rather than the value 
e k, l.l/// r reported bylPapaloizou & LarwoodI (120001) . In agree- 
ment with ICresswell et all ( l2007h . we find the peak positive an- 
gular momentum exchange rate (from the disc to the planet) to 
occur when e ~ 4H/ r, ra ther th an the value e ~ 2H/ r reported by 
IPapaloizou & LarwoodI (l2000l) . This may be an effect of adopt- 
ing a different surface density profile. 

Figure |2| shows the orbital evolution of a planet on a variety 
of inclined orbits, and the accompanying migration. As with the 
eccentric orbits, agreement between the hydrodynamic model 
and the modified A^-body code is generally good, and poorest 
at the transition to exponential decay. 

Figure [3] shows the orbital evolution of a planet on a variety 
of orbits which are both eccentric and inclined to the disc mid- 
plane. Agreement between the two numerical schemes is again 
good overall. The oscillations present in the most strongly ex- 
cited runs are physical artefacts, and possess a frequency re- 
lated to the precession rate of the longitude of ascending node. 
A bending distortion of the disc by the protoplanet is probably 
responsible, but does not significantly influence the underlying 
damping rate. It should be noted that our prescription for incli- 
nation evolution given by Eq.[T2|does not include this oscillatory 
behaviour. 

We note that it is not feasible to test planetary orbits with 
very high eccentricity or inclination against hydrodynamic mod- 
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Fig. 2. {Top) Inclination evolution of a 10 protoplanet on 
a variety of inclined orbits. Solid lines are the results from 
NIRVANA, crosses from the A^-body code. {Bottom) The cor- 
responding migration rates. 
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Fig. 3. {Top) Eccentricity evolution of a 10 Me protoplanet on 
a variety of eccentric and inclined orbits. Solid lines are the re- 
sults from NIRVANA, crosses from the A^-body code. {Bottom) 
Inclination evolution of the same orbits. 



els, due to the large range of scales over which the disc model 
must be simulated. However, previous experiments (Cresswell 
12006 ) and Paper I imply that swarms of protoplanets will spend 
the majority of their time with fairly low eccentricities and in- 
chnations. Equations [TTHlJ] successfully capture the character 



of the dynamical evolution {dy/dt oc f at low y and oc y^-^ at 
y» Hjr, for y = e or /) during the brief times when these quan- 
tities become large. 



4. Initial conditions and units 

We adopt a similar strategy to Paper I when setting up plane- 
tary initial conditions. We first define the semi-major axis of 
the innermost body to be a i - 5 AU, such that the popula- 
tion exists beyond the snow line where the proportion of solids 
available for protoplanet core formation is higher Successive 
protoplanet semi-major axes were determined by choosing their 
separations to be a specified number of mutual Hill radii. Thus 
- Qj + A^mH^mH whcrc A^mH IS typically 5 (though other val- 
ues are considered). The mutual Hill radius is defined by 



^mH — 



m, + mj\^ /flj + a 
3M. 



(17) 



For two planets on initially circular orbits, rapid instability oc- 
curs if the separation A between planets is less than the critical 
value 



^mH 



3.46 



(18) 



(iGladmanI [1993; Chambers et a l.l Il996h. Simulations of oli- 
garchic growth ( K okubo & Ida«200oF Thommes et al.l2003l) sug- 
gest that the mutual separation of protoplanets is normally ^ 
8 RmH- We typically adopt smaller spacings to maximise close 
encounters, although larger separations are also investigated. 

Planetary eccentricities were determined by defining a mean 
eccentricity jj^ for the planetary swarm, and standard deviation 
(Te - 0.01, with eccentricities then chosen randomly according 
to a Gaussian distribution. Each body is given a random argu- 
ment of pericentre. Inclinations are likewise Gaussian distributed 
according a mean yu, - 0° and standard deviation cr,, with ran- 
dom longitude of ascending node. 

Two different schemes were used to determine initial plan- 
etary masses. In the 'ordered' models, the standard procedure 
was to define the mass of the innermost body (usually m\ - 2 
Me) with subsequent bodies having m,+i = m, + 2 M^. This 
somewhat artificial set-up was chosen to maximise convergent 
migration, and hence maximise interactions between the bod- 
ies. In the 'randomised' models, a more natural distribution is 
formed from a randomly selected Gaussian distribution of plan- 
etary masses, with mean /i,„ and standard deviation cr,,,, subject 
to a lower-mass cutoff of either 2 or 0.2 M®, and an upper cutoff" 
of 20 Me (but note that collisions may raise masses above this 
value). 

The disc was initialised with scale height H/r - 0.05 and 
S(r) - So'" where So was typically chosen such that the 
disc contains 40 Jupiter masses of gaseous material within 40 
AU of the central star Other disc masses and surface density 
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Table 1. A subset of the 3D models performed. From left to 
right, the columns give: class name; mean eccentricity; standard 
deviation of inclination; mean mass; standard deviation of mass; 
lower-mass cutoff; initial separations in mutual Hill radii; disc 
mass (normalised against fiducial value); number of protoplan- 
ets. 



profiles were also used, and are described where appropriate 
in subsequent sections of this paper In the NIRVANA simula- 
tion the disc had a viscou s alpha parameter of or = 5 x 10"^ 
jShakura & Sunvaevll 19731) . 

The distribution of semi-major axes, planet and disc masses, 
and values of yu^, cTj. and cr, together define a class of model. For 
each model class five realisations of the initial data were gen- 
erated by rotating the random number seeds, giving rise to five 
different simulations. In total over 300 simulations were run; de- 
tails of a subset of the models used are given in table[Tl including 
those described in subsequent sections. 

5. Results of the hydrodynamic simulation 

One 3D hydrodynamic run (case HI) was performed to comple- 
ment the A^-body models. We adopted the fiducial ordered mass 
set-up, to maximise convergent differential migration. To model 
the system for over 2x10"* years while preserving the grid resolu- 
tion, eight protoplanets were used so that the radial extent of the 
grid could be lessened (ten protoplanets were normally used in 
the A^-body simulations described later). The protoplanets take 
the same masses as the eight inner bodies of the fiducial ordered 
runs, so that the inner body is 2 and the outer is 16 M^. The 
number of active cells used was (Nr,Ne,N^) - (288,40,444); 
this provides a resolution finer than those used to construct 
Eqns. fTTHlJj bv a factor of 1.5 in azimuth and 2 in radius. As in 
those tests, for the potential softening we adopt a value almost 
one tenth the vertical height of a cell. 

The results of this model are shown in Fig.lD Resonant mi- 
gration dominates the model: after a short period of low-level 
scattering and orbital exchanges (where horseshoe interactions 
rearrange the radial ordering of two adjacent protoplanets), the 
system settled down into two groups of planets which are in mu- 




5.0x10= 1,0'10° 1.5-10' 2.0-10° 
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Fig. 4. Orbital elements for the 3D NIRVANA model, similar to 
Fig. |5]but with eight protoplanets (run HI). Behaviour is sim- 
ilar to the early stages of the A^-body models, with scattering 
and collisions among the inner population. A co-orbital system 
is formed, which remains stable for the remainder of the inte- 
grations. Note only the smallest protoplanet achieves e > 0.1 
or / > 0.5° after the initial values have been damped. The 6 & 
10 Me protoplanets collide at f = 6.7 x lO"' yrs, and the 4 & 8 
Me protoplanets collide at f = 1.20 x 10"* yrs. 



tual mean-motion resonances, with all resonances being either 
4:3 or 5:4. The inner pair have suffered collisions (at 6.7 x lO-' 
and 1.2 x lO** yrs) and end the scattering phase in a 5:4 MMR 
and on orbits diverging from the external protoplanets due to 
their increased mass; they are no longer plotted once the inner 
body of the pair enters the inner damping region, but based on 
all other models performed, we expect this pair would continue 
to migrate inward in resonance. 

Resonant migration is thus the dominant outcome of the sim- 
ulation. Only the smallest, scattered body achieves a significant 
eccentricity or inclination, with the other protoplanets limited to 
e < 0.1 and / < 0.5° once their initial values have been damped 
by the disc. By the end of the run all other protoplanet inclina- 
tions are near zero due to the strong damping. In the absence of 
a magnetospheric cavity, or other halting mechanism, we expect 
all protoplanets to migrate into the central star 

One further point of note is that a co-orbital pair of plan- 
ets are observed. Initially this involves the planets undergoing 
mutual horseshoe motions, but it is expected from the results of 
Paper I that the damping action of the disc will reduce the libra- 
tion width until tadpole orbits result. The co-orbital pair remains 
stable in simultaneous resonant migration with other bodies on 
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both interior and exterior orbits for the remainder of the integra- 
tion. 

Despite a total run time of over 1 .5 x 10'' cpu hours on a par- 
allel facility, due to the high cost of performing hydrodynamic 
simulations in three dimensions we are unable to continue this 
simulation for more than a few x 1 0"^ years, nor repeat other simi- 
lar simulations many times in order to gain a statistical overview 
of these chaotic systems. Instead, we perform numerous A^-body 
simulations using the HENC-3D code to examine these issues. 

6. Results of the modified A^-body simulations 

We have performed more than 300 modified A^-body simulations 
to examine the evolution of clusters of low-mass protoplanets 
embedded in a 3D protoplanetary disc, varying planet masses, 
orbital parameters, and disc masses. Despite all this variation 
in initial conditions, a number of simple trends are observed. 
We present a select few results to illustrate these trends. Further 
detail of the 2D analogues of these trends may be found in Paper 
I. 

6.1. A fiducial ordered mass N-body simulation 

We choose one particular class of model (class 01 in table [TJ 
to act as the fiducial case, against which other models are com- 
pared. As shown in Fig. |5]the model follows a typical pattern: 
differential migration initially causes the planets to move closer 
to each other In the outer half of the swarm, orbital exchanges, 
where two planets appear to swap semi-major axes in a horse- 
shoe motion, may occur. A similar situation is seen among the 
inner half; however, due to the larger mass ratios ot,+i/ot; be- 
tween neighbouring inner bodies such an exchange leaves the 
less massive planet with a slightly larger semi-major axis than 
the body it has just replaced, according to the simple ratio 
Aaj/Aoi ^ mj/nii. The smaller body is now closer to the next 
external mass than the previous planet was, making further such 
exchanges (themselves with larger m^/m,) more likely. In this 
manner the smallest masses in the population are often passed 
outward from planet to planet, gaining in eccentricity and in- 
clination in the process. This process typically culminates with: 
(i) collision with another body; (ii) co-orbital capture; (Hi) ejec- 
tion beyond the outer edge of the swarm. A rarer fourth out- 
come is capture in the Hill sphere of another body, forming a 
planet-moon or binary planet system. HENC-3D is not equipped 
to model such encounters accurately, which always result in a 
collision after a few xlO"^ years. The process may repeat for sev- 
eral of the smallest bodies, such that the initially innermost plan- 
ets finally constitute (in some ordering) the outermost planets of 
the swarm, with the original outermost bodies now leading the 
inward migration having 'pushed through' the swarm. 

As the population of planets converges due to differential 
migration, a series of mean-motion resonances (MMRs) forms 
between the planets (or subsets of them), forcing them to mi- 
grate inward in lockstep. Occasional bursts of dynamical insta- 
bility occur, leading to the removal of bodies by collision, co- 
orbital capture or ejection to the outer edge of the swarm, and 
this helps stabilise the population. The initial phase of scatter- 
ing smaller bodies ends within a few thousand years, and the 
disc rapidly damps eccentricies and inclinations developed dur- 
ing this active phase. During the ensuing resonant migration, as 
separations between the planets decrease, a single planet is oc- 
casionally perturbed sufficiently to break the resonant chain, and 
the process of the swarm reordering itself may repeat on a shorter 
time scale (often < lO-' yrs), as seen in Fig.|5]at time t x 6.5 x 10^ 
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5,0x10" 1.0x10^ 1,5x10^ 2,0x10^ 

time (years) 



Fig. 5. Evolution of a 10 planet A^-body model with ordered 
masses and the fiducial set-up (Ol). (Top) Semi-major axes 
of the migrating embryos. Short periods of activity are fol- 
lowed by long periods of migration between bodies in first or- 
der mean-motion resonances. The 4 & 18 M® planets collide at 
t - 6.76 X 10'* yrs. (Middle) The embryos' eccentricities over the 
same time. (Bottom) The embryos' inclinations. 



years. Ultimately the population, now comprised of an individ- 
ual planet and two resonant groups, migrates towards the cen- 
tral star, where it will be accreted in the absence of a stopping 
mechanism, such as a disc edge created by a magnetospheric 
cavity, or a 'planet trap' associated wi th a region of the d isc with 
a positive surface density gradient ( Masset et al.ll2006l) . The fi- 
nal masses of the planets, labelled from smaller to larger radii at 
the end of the simulation, are (first group) 8,22, 14,20, 16 
(second group, including the co-orbital pair) 10, 6, 12 Me, and 2 
Me, showing how the population has separated with the smallest 
bodies at larger radii and the largest bodies closer to the star. The 
smallest planets, if scattered significantly far outwards, or with 
sufficient eccentricity or inclination to slow their migration, may 
survive for times on the order of 10^ years, raising the possibil- 
ity of survival as a super-terrestrial body if the population formed 
late in the disc's lifetime. 

Due to the strong inclination damping, the system is approx- 
imately planar for most of its duration. As the planets migrate 
together, the MMRs between neighbouring bodies are almost al- 
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Fig. 6. The first order mean-motion resonances between the planets in Fig.|5] during the longest resonant migration phase. Each 
plot is labelled (Pi,Pj)-p:q, where P, and Pj are the exterior and interior protoplanets in each resonance respectively, and p:q is the 
commensurability between them. 



ways (> 99%) first order, typically 3:2, 4:3,. . . ,8:7. The first or- 
der and co-orbital resonances between the planets in Fig. |5] are 
shown in Fig.|6] where the planets are labelled 1-10 with 1 be- 
ing the initially innermost body, and 10 the initially outermost. 
The resonant angles for the first order resonances (p = q+l) are 
defined by: 

4>i - pAi - qA2 - Ci>i 

(p2 = pAi - qA2 - (jJ2 (19) 

where A\ (A2) and wi (W2) are the mean longitude and lon- 
gitude of pericentre for the outer (inner) planet, respectively. 
Commonly, the complicated lattice of resonances between 
groups of planets generates additional first, second and third or- 
der resonances between non-neighbouring planets {e.g. alternat- 
ing 5:4 and 6:5 resonances implies the existence of a 3:2 reso- 
nance between the non-neighbouring planets). Rather than being 



a source of perturbation to resonant pairs, additional planets thus 
often provide further stability to a migrating resonant group. 

In Sect. |6.5| we present a more statistical analysis of our re- 
sults, but in the next subsection we discuss how the qualitative 
nature of the results change as initial conditions are varied. 



6.2. Variation of tine initiai conditions 

The behaviour discussed above represents the typical behaviour 
observed among the ordered mass models. Many other sets of 
initial conditions were used, including varying initial eccentric- 
ities and inclinations up to p^ - 0.15 and cr,- - 6° (e.g. runs 
02-04). 
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6.2.1. Variation of initial eccentricities 

Increasing the initial eccentricities typically produces a much 
longer dynamically active phase among the swarm, which in- 
volves more planets and prevents the formation of early but tem- 
porary resonant groups. This is especially true when the initial 
eccentricities are such that protoplanets begin on crossing orbits. 
Much of the population is involved in dynamical relaxation until 
several bodies have been expelled from the group, or lost in col- 
lisions. In the long run, however, the final states of these systems 
are similar, and consist of inwardly migrating groups of planets, 
in multiple mean-motion resonances, often including co-orbital 
systems. Their final fate remains accretion by the central star. 

6.2.2. Variation of initial inclinations 

Raising initial inclinations has little effect, as they are damped 
by the disc before large-scale scattering activity begins. This is 
to be expected as raising inclinations does not lead to the cre- 
ation of crossing orbits, unlike the initial eccentricity variations. 
Systems with large values of initial inclination have their incli- 
nations damped on the time scale of a few thousand orbits, after 
which they enter a phase of dynamical interaction similar to that 
shown by the fiducial run described in Sect. 16.11 

Scattering throughout the population is now more common 
overall than in 2D (see Sect. 16. 5t and the resulting activity 
largely removes memory of the initial conditions. For the initial 
eccentricity/inclination distribution to substantially alter the long 
term evolution requires seemingly unphysical values for these 
quantities, corresponding to non-disc-like initial conditions. 

6.2.3. Variation of initial separations 

The initial separations between the planets were varied, with val- 
ues ranging between 4 and 8 ^mH- This was found to have no 
qualitative effect on the final states of the systems. Early stages 
showed significant differences, with closer protoplanets under- 
going prolonged periods of violent excitation, and more distant 
ones slowly moving closer together and taking longer to pro- 
duce any significant scattering activity. Once scattering activity 
ceased to continue the final states of the systems were again very 
similar to those already described. 

6.2.4. Modifying the disc mass 

Altering the disc mass produced more substantial changes. 
Reduction by a factor of 2 (class 07) produced population-wide 
gravitational interactions for similar lengths of time (~ 1-2x10^ 
years) despite longer migration times, yet the swarm was more 
likely to retain its original order, with fewer of the more massive 
planets passing through the population and instead remaining at 
the rear of the population where they drive the whole group for- 
ward. This presumably arises because of the overall reduction 
in convergent migration rates. The resonant groups were often 
larger, comparable to those seen in 2D (consisting of up to 7 
bodies). However, those small planets that were scattered tended 
to achieve higher eccentricities (e ^ 0.5-0.6) and higher incli- 
nations (6°-12°), and remained in such excited orbits for several 
times longer than in the fiducial case. By spending significant 
portions of each of orbit away from the other protoplanets (due 
to the high inclinations), encounters that may result in collisions 
or co-orbital pairs become less likely, allowing the smaller body 
to retain its excited state for longer Further reductions in disc 



mass increased the effect slightly and were more likely to pro- 
duce larger resonant groups. 

6.2.5. Changing the number of planets 

We have run models with smaller numbers of protoplanets (five 
rather than ten). One such model is listed as class 09 in table[T] 
Although the early phase of evolution can occasionally involve 
all planets being in mutual mean-motion resonances, this config- 
uration was found to be stable for long time scales in only ap- 
prox. 20% of runs; interestingly, only those models that formed a 
stable co-orbital pair were able to sustain such a five-planet res- 
onant group. In all other cases instabilities leading to scattering 
reduced the group size, either by ejecting one of the smallest pro- 
toplanets to a larger orbit (most common when the initial mass 
range of protoplanets was largest, 2-20 M®) or by removing one 
or more protoplanets by collision (most common when the initial 
mass range was smallest, 2-10 Me). In all cases the long-term 
evolution always resulted in migration of the remaining proto- 
planetary swarm into the central star. 

6.3. A fiducial randomised mass N-body simulation 

When initiating models with randomised masses, we set a mean 
fi,„ and standard deviation cr,„ for the mass distribution, and 
chose masses randomly according to a Gaussian distribution. A 
wide range of mass distributions were considered, which were 
conflated with the variations in initial conditions used in the or- 
dered mass models. For a fiducial randomised model, we choose 
initial conditions and the planetary mass range similar to the 
fiducial ordered case (class Rl). Such a model is shown in Fig.|7] 

Overall, randomised models follow broadly similar evolu- 
tionary paths to the ordered models, except that differential mi- 
gration carries some planets, or groups of planets, away from 
each other. Consequently, unless the early gravitational interac- 
tions among the population are prolonged or involve especially 
strong scattering, the protoplanets rapidly (< 5 x 10^ yrs) break 
up into smaller groups of typically 2-5 bodies, with the slowest 
migrators typically following behind in isolation after scatter- 
ing. This behaviour is clearly seen in Fig. [T] and also in Fig. [8] 
which shows snapshots of the end-states of a random selection 
of models. Planets of the largest or smallest mass typically end 
the simulation migrating alone or in small groups at the inner or 
outer edges of the population, respectively. 

6.3.1 . Variation of the mass distribution 

Several different mass distributions were combined with the var- 
ious initial eccentricity and inclination distributions that we have 
considered previously. Values in the range 5 Me < < 10 
Me and 3 Me < cr,,, < 7 M® (e.g. classes R1-R5) were used in 
setting up initial protoplanetary swarms. The resulting models 
can be characterised by a number of features. 

First it was found that among the primary mass ranges stud- 
ied (5 < /im < 10 Me), when cr,„ ^ /^,„/3 scattering was substan- 
tially reduced, resulting typically in one large group migrating 
without further strong interaction after an initial period of or- 
bital readjustment. This is clearly due to the fact that the planets 
are migrating inward at similar rates. For larger ratios of cr„,/yu,„ 
planet-planet scattering and its natural consequences (collisions, 
co-orbital planets, etc.) were more common, due to the concomi- 
tant increase in differential migration. Once strong scattering 
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Fig. 7. Evolution of a 10 planet A^-body model with randomised 
masses and the fiducial set-up (Rl). (Top) Semi-major axes of 
the migrating embryos. Short periods of activity are followed by 
long periods of migration of bodies in first order mean-motion 
resonances. (Middle) The embryos' eccentricities over the same 
time. (Bottom) The embryos' inclinations. 



was induced, however, we found no strong correlations between 
the final outcomes and the value of cr„,/yu,„. 

The second feature we note is that additional tests for lower- 
mass planets (see Sect. 16.4.11 for a more complete discussion), 
using mean masses in the interval (0.3 ^ ;U,„ ^ 2 Me), seem to 
result in more energetic and population-wide scattering activity 
than is the case for higher-mass planets, which may seem sur- 
prising at first glance. This result appears to originate in the fact 
that the lower-mass planets tend to become trapped in first order 
resonances of high degree (even as large as 14:13 for the lowest 
mass cases considered), placing the planets in very close prox- 
imity. Localised instabilites in the resonantly migrating swarm 
then lead to strong dynamical interaction and scattering. Such 
large-scale activity was also seen among ordered mass runs us- 
ing similarly small initial masses. 

The third feature shown by the randomised mass runs is that 
the initial distribution in semi-major axes of the bodies deter- 
mines much of the subsequent evolution. If the innermost bodies 
are of a higher mass than the mean of the population, those plan- 
ets will form a separate group and migrate away from the rest. 
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Fig. 8. Snapshots of two randomly selected simulations from the 
model classes shown, with the exception of the first (Rl) and 
fourth (R2) rows, which correspond to Figs.|2]&|9]respectively. 
The snapshots are taken at f = 2 x 10^ yrs. The area of each 
planet is linearly proportional to its mass; the smallest body (in 
the top row) is of 2 M®, and the largest (formed from multiple 
collisions) is of 39 M®. The co-orbitals of Fig.|9]are clearly vis- 
ible as concentric circles. 



generally leading to an uneventful inward journey for the whole 
population. Conversely, if several massive bodies are located in 
the outer half of the population, then they are likely to drive a 
larger resonant group ahead of themselves, usually causing the 
smallest bodies to be scattered outwards, and also more likely to 
generate co-orbital systems. 

The fourth feature we observe is that the initial distribution 
of masses has more effect on the resulting level of activity than 
their initial separations. That is, although it may take longer for 
differential migration to bring together initially widely separated 
bodies, the resulting activity is determined by the mass distribu- 
tion just as it is when the protoplanets are initially closer to- 
gether Separations of up to 20 RmH have been tested, with the 
mass distribution remaining the dominant factor 

Together, these features suggests that the level of activity 
among a population of moderate to large cores, assuming a mod- 
erate spacing of several mutual Hill radii, is dependent on the 
mass ratios of nearby bodies, with little scattering activity be- 
low a critical value cr„,//i„, $1/3 for a Gaussian distribution. 
Collections of smaller masses appear to undergo periods of evo- 
lution where the degree of scattering activity is greater than for 
their more massive counterparts. We note that the smaller pro- 
toplanets generally form MMRs of higher degree, in the range 
8:7-10:9 (or in more extreme cases 14:13), which in conjunc- 
tion with the weaker damping forces on eccentricity and incli- 
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nation may account for the increased scattering due to closer 
proximity. Nonetheless, we observed that in some cases these 
resonances could remain stable over long periods, and examples 
of long term stable planets in the 14:13 resonance occurred in 
some simulations. 

Varying the mass distribution produced no significant corre- 
lation with any other parameter variation (Sect. | 



6.4. Prevalence of co-orbital systems 

A significant feature that arose in Paper I was the unexpected 
abundance of stable co-orbital planets, orbiting around their mu- 
tual L4/L5 points. Co-orbital systems form when planet-planet 
scattering causes a planet to be perturbed such that its semi- 
major axis becomes very similar to that of another planet in the 
system. Disc-induced eccentricity damping then ensures rapid 
decay of the planet eccentricity, leading to co-orbital capture. 
Within the planetary swarm, the condition for long-term capture 
is simply that this eccentricity decay occurs before the scattered 
planet undergoes a close encounter with another protoplanet, 
which would otherwise disrupt the co-orbital system. We find 
that many of the co-orbital systems remain stable for the du- 
ration of the simulation while migrating inward over distances 
greater than 10 AU. 

In 3D, we find that these co-orbital planets are more common 
than in 2D (see Fig.fTOli. occurring in almost 45% of ordered and 
almost 35% of randomised simulations; of these, approximately 
21% and 14% respectively contained multiple examples of sta- 
ble, co-orbital pairs in the same simulation (as with other res- 
onances, co-orbital pairs are deemed stable within a simulation 
if they survive for a time > 10^ years.) We note that the ec- 
centricity damping rate adopted in the 2D simulations of Paper 
I was smaller than that used for the 3D runs in this paper. We 
attribute this increased co-orbital frequency to the higher ec- 
centricity damping rate and the generally smaller sizes of the 
resonant groups — both of which reduce the opportunity for a 
recently captured co-orbital body to be disturbed by other pro- 
toplanets — together with a slight increase in overall scatter- 
ing activity (Sect. [63] ). Inclined orbits are not an obstacle to co- 
orbital formation, and planets with inclinations greater than 10° 
are readily captured into stable horseshoe orbits, with these mu- 
tual inclinations eventually being damped by the disc. 

When a co-orbital system first forms it is usually in a mu- 
tual horsehoe orbit. In all except one instance the horseshoe mo- 
tions decayed because of the disc's action to form tadpole orbits, 
maintaining small oscillations about the L4/L5 points. The transi- 
tion from horseshoe to tadpole motion typically takes 0.3-1x10"* 
years for our standard disc parameters. 

A detail of two migrating co-orbital systems (from class R2) 
is shown in Fig.|9l Initially the horseshoe librations are of large 
amplitude, but these decrease with time, and after 2-3 xlO"' yrs 
tadpole motions result. The system then contines to migrate in- 
ward maintaining this architecture. 

Co-orbital planets are found both as isolated pairs, and as 
part of larger resonant groups, sharing mean-motion resonances 
with external and internal bodies like any other individual body. 
At the end-state of a simulation, a co-orbital pair may thus be 
in multiple first- and second-order MMRs with two or more in- 
terior/exterior bodies, all densely packed within 0.4 AU of the 
star. 



1 1 1 1 1 

t = 74873 yrs - 
• 


"'"--A... 


1 1 1 


A \ 




m 











-4 


-2 


2 


4 


t = 99886 y 


»•-■ ' 

a ' 


• 










A 


, 1 


1 








-4 -? 



■ t = 149B50 yrs 


















. . 1 ... 1 


. . . 1 . 





Fig. 9. Evolution of a co-orbital system (class R2) migrating in- 
wards: two co-orbital pairs, locked in a 3:2 MMR. The size of 
each planet is proportional to its mass, while the open triangles 
display each primary body's L4 & L5 points. The other proto- 
planets present in the model are not displayed, but those shown 
are part of a resonant group encompassing nine bodies. (Top) 
Initially the co-orbital planets' motions cover the whole horse- 
shoe region. (Middle) The disc's action causes the planets to shift 
into tadpole orbits. (Bottom) The planets migrate inwards under 
small librations. 



6.4.1. Limits on co-orbital formation 

Other studies have examined the behaviour of swarms of pro- 
toplanets, subject to type I migration, at earlier stages of for- 
mati on when protopla netary masses are considerably smaller 
(e.g. iMcNeil et al.ll2005i) . Similar patterns of resonant migration 
have been observed, but co-orbital planets were not found, in 
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contrast to their ubiquity here and in Paper I. Consequently, we 
seek to determine limits on co-orbital formation and explain this 
discrepancy between models; we focus on planetary mass, initial 
separations, and the surface density profile of the disc. 

To compare with studies conducted using lower-mass proto- 
planetary swarms, we selected one fiducial class from each of the 
ordered (Ol) and randomised (Rl) models, and reran them with 
the protoplanets' masses (yu,„ and cr,„, in the randomised case) 
repeatedly halved; the lower-mass cut-offs are set to in both 
cases (see Sect. 16.3.11 ). Additionally, in the randomised models 
we increased the initial separations in stages up to 20 RmH- 

When using our standard disc model with these smaller pro- 
toplanet masses, 25% of simulations resulted in stable co-orbital 
pairs for both the ordered and randomised mass distributions, 
averaged over the ranges of planetary masses and initial separa- 
tions considered. 

With smaller sample sizes, we can be less confident about 
these statistics than for the main body of the simulations per- 
formed. However, the fact that co-orbitals continue to appear at 
all, and with frequencies of ~ 25%, is significant. The region 
within which co-orbital capture can occur scales ~ nip^^, and so 
we may expect capture rates to fall to near zero for populations 
of such small bodies. However, differential migration brings pro- 
toplanets closer together; as masses decrease, the MMRs formed 
between neighbouring protoplanets are of higher degree (sta- 
ble resonances up to 14:13 were observed among this work, 
as mentioned previously) and so the planets lie closer together. 
Consequently only a small shift in the orbit is necessary for one 
body to enter the horseshoe region of another and be captured, 
and a perturbation from a third body can provide the necessary 
energy to jump into the horseshoe region. All that is required for 
co-orbital formation is 3-4 protoplanets in close proximity. 

We also note that, for the randomised mass distributions, co- 
orbital frequency is essentially independent of the initial sep- 
arations. This is easily interpreted as resonant migration driv- 
ing groups of bodies together; differential migration may cause 
a relatively massive protoplanet to 'sweep up' several smaller 
bodies in MMRs ahead of it. Formation of a co-orbital pair 
from such collections is then dependent only on the usual mu- 
tual interactions within a small group. Provided the initial sep- 
arations are not so large that such groups do not have time to 
for m — highly un l ikely for most oligarchic growth scenarios, 
e.g. Kokubo & Idal ( 1200 0) — co-orbital formation thus remains 
primarily dependent on the initial mass distribution^ 

On e feature different here from the set-up of iMcNeil et al.l 
(l2005 l) are the disc damping times, where they utiU sed the higher 
damping rates of Papaloiz ou & LarwoodI (l2000l) ; the discrep- 
ancy is approximately given by a factor of 3 (see also Figs. 1 
& 2 of Paper 1). To test the effect of this difference, we shortened 
the damping time in Eqn. [TT]by this factor, and reran a selection 
of models with the lowest mass distributions previously consid- 
ered. We find that although interactions between the planets are 
significantly reduced, co-orbital planets are still able to form as 
differential migration brings a population together, with crowd- 
ing producing the minor perturbations and individual orbital ex- 
changes required, despite low eccentricities among the group. 

The situation changes dramatically when the disc surface 
density profile is steepened. We ran a suite of models with the 
power law exponent being decreased to -3/2, which is the value 
usual ly adopted for the minimum mass solar nebula jHavashil 
Il981h . and corresponds to the disc model used in McNeil et al 
(2005). In this case the incidence of co-orbital system forma- 
tion fell to almost zero, this being a direct consequence of the 
fact that the migration of interior bodies is speeded up, and that 



of bodies lying further out in the disc is slowed down, when 
the density profile is steepened. Thus, it would appear that co- 
orbital planet formation is highly sensitive to there being strong 
differential migration that brings bodies together, a situation that 
is highly favoured in discs with flatter surface density profiles. 
Observation of co-orbital extrasolar planets would thus be an 
indicator of strong dynamical relaxation occuring during planet 
formation, induced and maintained in-part by there having been 
a suitably flat radial surface density profile in the original proto- 
planetary disc. 

6.4.2. Long-term stability 

To test the stability of these tadpole planets during and beyond 
disc dispersal, a time-dependent function was added to Eq.|9]that 
reduced the disc-induced forces acting on the planets exponen- 
tially, simulating disc mass loss according to a prescribed time 
scale. We then selected several simulations that produced tad- 
pole planets, removed the other planets, and restarted the model 
after the formation of the co-orbital pair. Each pair was found to 
separate if the disc forces were reduced too rapidly, typically cor- 
responding to a (grossly unphysical) mass-halving time of under 
lO^* years. In all other cases, the co-orbital pair remained stable, 
with their orbits remaining largely unchanged once the remain- 
ing simulated disc mass became negligible. 

Test calculations were also run in which the non-co-orbital 
planets were not removed from the system, and these contained 
instances where: the co-orbital pair was largely isolated from 
all other bodies because of prior differential migration; the co- 
orbital pair were on significantly eccentric orbits (both planets 
possessing e > 0.25); the co-orbital pair were in resonance 
with additional bodies driving faster inward migration (while 
the migration force remained effective). Subject to the condition 
on disc mass-halving time stated above, all these models were 
found to be stable for as long as the integrations were continued, 
with a minimum simulation time of ~ 2 x 10^ years in each case 
and a typical duration an order of magnitude greater. One model 
was allowed to evolve in a disc-free environment for over 4x10^ 
years, with the co-orbital system remaining stable for this time. 

We note those models with additional planets in mean- 
motion resonance with the co-orbital pair also implies the long 
term stability of pairs of protoplanets in MMR after disc disper- 
sal, and these resonant systems often contain numerous bodies in 
mutual MMRs. While the observed exoplanets in MMR are gi- 
ant planets (e.g. GJ 876, 55 Cancri, HD 12831 1, HD 82943, HD 
7352), whose resonances are thought to have been established 
through differential type II migration (e.g. ISnellgrove et al.l 
l2001l;lLee & Pealel200alKlev et al.l2004l) . our simulations show 
that differential type I migration may lead to the formation of 
multiple MMRs between groups of lower-mass planets. These 
will become amenable to detection as observational techniques 
allow greater exploration of the lower-mass end of the extrasolar 
planet population. 

6.5. Statistics ofsimuiation outcomes 

We now discuss the frequency with which certain simulation 
outcomes, such as planet-planet collisions, co-orbital system for- 
mation, etc. arose. Figure[TO]displays the frequency with which 
co-orbital systems formed within the simulations, and survived 
for at least 10^ years. The two leftmost bars show results for 
the 2D simulations presented in Paper I. The rightmost bars 
show results for the 3D runs presented in this paper. In the 
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Fig. 10. Frequency of A^-body models displaying stable horse- 
shoe or tadpole planets. On the left, A^-body data from Paper I is 
shown for comparison. Left-right: 2D ordered, 2D randomised, 
3D ordered, 3D randomised. 




Fig. 11. Frequency of collisions per simulation for the ordered 
mass models (left) and randomised mass models (right). The 
open bars show the corresponding 2D data of Paper I. 



cases where the protoplanet mass increased with increasing ini- 
tial semi-major axis, we see that the 2D and 3D results are very 
similar, with ~ 40% of each set of calculations producing long- 
lived co-orbital systems. The situation for randomised mass dis- 
tributions is different, however, with only 3% of 2D runs lead- 
ing to co-orbital formation, whereas 34% of 3D runs resulted 
in long-term co-orbital formation We believe this is due to a 
stronger eccentricity damping rate (in line with lTanaka & WardI 
|2004|) in the 3D models, which allows the orbits of recently cap- 
tured co-orbital bodies to be circularised and achieve a tadpole 
orbit sooner, reducing the probability of the body encountering 
a third protoplanet and being scattered from its horeshoe orbit. 
The effect is lessened among the models with initially ordered 
mass distributions because these favour large groups of proto- 
planets, meaning such a scattered body is likely to encounter 
further protoplanets and have multiple opportunities to be cap- 
tured as a co-orbital entity. 

Figure [TT] shows the collision frequency for 2D and 3D sim- 
ulations, as a function of the number of collisions occurring per 
simulation. The left hand charts show results for runs where 
the initial protoplanet mass increased with semi-major axis, the 
right hand charts show results for randomised mass distributions. 
Although there are small differences between both 2D and 3D 
runs, and between ordered and randomised mass distributions, 
the striking feature of these plots is their similarity. One nor- 
mally expects that 2D and 3D simulations would produce dif- 
ferent collision frequencies (with 2D runs leading to many more 
collisions), but this is not borne out by our results. The reason is 
most likely to be that the very strong inclincation damping pro- 
vided by the disc causes the planetary swarms to remain quasi- 
2D, thus increasing the collision frequency. 

Figure [12] shows the frequency with which resonant groups 
form as a function of the number of protoplanets contained in 
the resonant group for 3D simulations only. Resonant groups are 
only counted in simulations if they survive for 10^ years or more. 
The left hand chart is for models in which the initial planet mass 



increased with semi-major axis, and that on the right hand side 
is for the randomised mass distributions. The 3D runs show that 
simulations containing resonant groups of just two planets are 
most common, but with a significant number containing three to 
six planets. This is a clear indication that resonances involving 
smaller numbers of planets are more stable over the longer term. 

7. Conclusions 

We have performed simulations, using two different numerical 
schemes, which examine the evolution of swarms of low-mass 
protoplanets embedded in a 3D protoplanetary disc. First, we 
used hydrodynamical simulations to model the orbital evolution 
of planets on initially eccentric and/or inclined orbits, and fit- 
ted analytic equations for the rate of eccentricity and inclina- 
tion damping, and migration, to these simulation results. These 
equations were then incorporated into an A^-body code, which 
was used to perform many simulations of planetary swarms em- 
bedded in protoplanetary discs. We also performed a single hy- 
drodynamic simulation of a planetary swarm embedded in a 
disc to verify the qualitative outcome of the A^-body simula- 
tions. Although this simulation could not be evolved for as long 
as the A^-body runs due to computational cost, we found basic 
agreement between it and the A^-body simulations for the first 
^ 30, 000 years of evolution. 

The main aim of this work was to re-examine previous re- 
sults we had obtained using 2D simulations (see Paper I). It 
is known that type I migration of low-mass protoplanets may 
be slowed or even reversed by sustaining significant eccentrici- 
ties (Papaloizou & Larwood 2000), and in Paper I we examined 
whether gravitational interactions between a swarm of proto- 
planets can provide the necessary excitation of eccentricitity to 
prevent inward migration. The conclusion of that work was that 
disc-induced eccentricity damping is too strong, and ultimately 
protoplanetary swarms migrate into the central star In this paper 
we have examined whether the inclusion of 3D effects changes 
this basic conclusion, as the dynamics in 3D can be quite differ- 
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Fig. 12. Frequency of resonant group size per simulation for the 
ordered mass models (left) and randomised mass models (right). 



ent. For example, the collision frequency may be reduced, and 
the possible excitation of significant inclination may assist in 
slowing down migration. 

Our results indicate that the collision rate is reduced only 
marginally in 3D, because the strong disc-induced inclination 
damping causes the system to remain quasi-2D. We find that the 
disc-induced eccentricity damping remains too strong, so that 
a protoplanetary swarm is unable to maintain long epochs of 
strong gravitational scattering with concomitant high eccentrici- 
ties. 

Despite a wide variety of initial conditions, a typical mode 
of behaviour is observed across the models, which may be sum- 
marised as follows: Differential type I migration leads to con- 
verging orbits and a crowded system. Gravitational scattering 
leads to orbital exchanges between neighbouring planets, for- 
mation of groups of planets in mutual mean-motion resonances, 
formation of 1:1 co-orbital resonances, and occasional colli- 
sions. Smaller protoplanets are frequently scattered out beyond 
the population, but rarely achieve a semi-major axis greater than 
that of the initially outermost protoplanet prior to migration, nor 
an eccentricity or inclination capable of significantly prolong- 
ing their life times. Within < 10^ yrs of commencement, the 
simulated system has typically settled down into a state with 
the protoplanets being in resonant groups, with each member of 
the group being in a mutual mean-motion resonance. These are 
usually first order mean-motion resonances, typically of degree 
3:2-8:7. These resonant groups migrate inward in lockstep, to be 
accreted by the central star in the absence of a stopping mech- 
anism such as an inner magnetospheric cavity. Occasionally a 
late burst of scattering activity occurs once a swarm of four or 
more bodies has migrated over several AU, but these systems 
always settle down to another phase of resonant migration that 
takes them to the star A small, slowly migrating protoplanet may 
sometimes 1% of runs) be scattered to the outer edge of the 
swarm (where it cannot be trapped in resonance by a faster mi- 
grating body) and survive for ~ lO*" yrs. 

We conclude that if multiple protoplanets form coevally from 
oligarchic growth in the giant planet zone of a laminar disc, then 



the long term evolution of the system will usually be collective 
inward migration and ultimate accretion by the central star This 
occurs on a time scale much shorter than the accretion time re- 
quired to accrete gaseous envelopes and end ty pe I migration 
jPoUack et al1ll996l:lPapaloizou & Nelsonll2005l) . 

Co-orbital planets form as a natural consequence of gravita- 
tional scattering in crowded systems, and occur in more than one 
third of all models we have considered. Simulations performed 
with sub-Earth mass protoplanets showed that even very low- 
mass bodies could form and maintain co-orbital systems. We 
showed that an important factor in establishing these systems 
was the existence of a relatively flat disc surface density profile, 
which promotes convergence of planetary orbits through differ- 
ential type I migration. We also showed that these co-orbital sys- 
tems can be stable for over 3 billion years after gas disc dispersal. 
Co-orbital triples (i.e. three planets all in mutual 1:1 resonance) 
are also potentially viable, though none were found in the suite 
of simulations presented in this paper. Previous A^-body simu- 
lations, using a slightly modified form for the eccentricity and 
damping formulae adopted in this paper, did result in such sys- 
tems occasionally forming (Cresswell 2006). Although no plan- 
etary system is known to host such a configuration, we note that 
such 'tadpole twins' have been observed in orbit around Saturn, 
with H elene and Polyduec es occupying the L4 and L5 points of 
Dione (Murra v et alj|2005l) . The observation of co-orbital extra- 
solar planet systems would be a strong indicator that a period of 
sustained dynamical relaxation had occured during the forma- 
tion of that system. 

There are a number of open questions regarding the co- 
orbital planets that form in our simulations. The main issue 
is whether they can remain stable if one or both planets un- 
dergo gas accretion to become giants, since such a system 
would be more amenable to detection among the extrasolar 
planet population. Another question regarding our overall re- 
sults is whether a combination of planet-planet interactions 
and stochastic mig ration, induced by turbulent density fluctu- 
ations in the disc jNelson & Papaloizoul 120041; iNelsonI l2005h . 
can help prevent large-scale migration of some planets within 
the swarm. In particular, we find that mean-motion resonances 
are very effective in shielding the protoplanets from close- 
encounters that lead to scattering, and stochastic torques may 
increase the fragility of these resonances, including the co- 
orbitals. Another point of interest is the effect of other halting 
mechanisms for type I migra tion, such as a disc edge cau sed 
by a magnetospheric cavity dTerquem & Papaloizoul 120071) or 
a 'planet trap' produc ed by a region of positive d ensity gra- 
dient in the disc (M asset et all l2006t iMorbidelU e t al. 2008|), 
on a resonant group that has already for med and migrated 
some distance. iTerquem & Papaloizoul (l2007i) suggest that near- 
commeasurabilities between prot oplanets will survive bey ond a 
disc edge, while experiments by iPierens & NelsonI ( 12008 ) sim- 
ilarly indicate that a disc edge formed by the action of a close 
stellar binary companion can prevent the further migration of a 
moderate (~ 5) swarm of protoplanets, possibly subject to fur- 
ther reordering of the population. Our results suggest that the 
survival of co-orbital systems in such circumstances is possible, 
but depends on the rate at which migration is halted and any 
further scattering among the system. We will address these and 
other issues relating to planetary swarms and co-orbitals in a fu- 
ture paper 
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